<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<meta http-equiv="X-UA-Compatible" content="IE=9"/>
<meta name="generator" content="Doxygen 1.9.1"/>
<meta name="viewport" content="width=device-width, initial-scale=1"/>
<title>Eigen: Sparse matrix manipulations</title>
<link href="tabs.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="jquery.js"></script>
<script type="text/javascript" src="dynsections.js"></script>
<link href="navtree.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="resize.js"></script>
<script type="text/javascript" src="navtreedata.js"></script>
<script type="text/javascript" src="navtree.js"></script>
<link href="search/search.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="search/searchdata.js"></script>
<script type="text/javascript" src="search/search.js"></script>
<script type="text/javascript">
/* @license magnet:?xt=urn:btih:cf05388f2679ee054f2beb29a391d25f4e673ac3&amp;dn=gpl-2.0.txt GPL-v2 */
  $(document).ready(function() { init_search(); });
/* @license-end */
</script>
<script type="text/x-mathjax-config">
  MathJax.Hub.Config({
    extensions: ["tex2jax.js", "TeX/AMSmath.js", "TeX/AMSsymbols.js"],
    jax: ["input/TeX","output/HTML-CSS"],
});
</script>
<script type="text/javascript" async="async" src="https://cdn.mathjax.org/mathjax/latest/MathJax.js"></script>
<link href="doxygen.css" rel="stylesheet" type="text/css" />
<link href="eigendoxy.css" rel="stylesheet" type="text/css">
<!--  -->
<script type="text/javascript" src="eigen_navtree_hacks.js"></script>
</head>
<body>
<div id="top"><!-- do not remove this div, it is closed by doxygen! -->
<div id="titlearea">
<table cellspacing="0" cellpadding="0">
 <tbody>
 <tr style="height: 56px;">
  <td id="projectlogo"><img alt="Logo" src="Eigen_Silly_Professor_64x64.png"/></td>
  <td id="projectalign" style="padding-left: 0.5em;">
   <div id="projectname"><a href="http://eigen.tuxfamily.org">Eigen</a>
   &#160;<span id="projectnumber">3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)</span>
   </div>
  </td>
   <td>        <div id="MSearchBox" class="MSearchBoxInactive">
        <span class="left">
          <img id="MSearchSelect" src="search/mag_sel.svg"
               onmouseover="return searchBox.OnSearchSelectShow()"
               onmouseout="return searchBox.OnSearchSelectHide()"
               alt=""/>
          <input type="text" id="MSearchField" value="Search" accesskey="S"
               onfocus="searchBox.OnSearchFieldFocus(true)" 
               onblur="searchBox.OnSearchFieldFocus(false)" 
               onkeyup="searchBox.OnSearchFieldChange(event)"/>
          </span><span class="right">
            <a id="MSearchClose" href="javascript:searchBox.CloseResultsWindow()"><img id="MSearchCloseImg" border="0" src="search/close.svg" alt=""/></a>
          </span>
        </div>
</td>
 </tr>
 </tbody>
</table>
</div>
<!-- end header part -->
<!-- Generated by Doxygen 1.9.1 -->
<script type="text/javascript">
/* @license magnet:?xt=urn:btih:cf05388f2679ee054f2beb29a391d25f4e673ac3&amp;dn=gpl-2.0.txt GPL-v2 */
var searchBox = new SearchBox("searchBox", "search",false,'Search','.html');
/* @license-end */
</script>
</div><!-- top -->
<div id="side-nav" class="ui-resizable side-nav-resizable">
  <div id="nav-tree">
    <div id="nav-tree-contents">
      <div id="nav-sync" class="sync"></div>
    </div>
  </div>
  <div id="splitbar" style="-moz-user-select:none;" 
       class="ui-resizable-handle">
  </div>
</div>
<script type="text/javascript">
/* @license magnet:?xt=urn:btih:cf05388f2679ee054f2beb29a391d25f4e673ac3&amp;dn=gpl-2.0.txt GPL-v2 */
$(document).ready(function(){initNavTree('group__TutorialSparse.html',''); initResizable(); });
/* @license-end */
</script>
<div id="doc-content">
<!-- window showing the filter options -->
<div id="MSearchSelectWindow"
     onmouseover="return searchBox.OnSearchSelectShow()"
     onmouseout="return searchBox.OnSearchSelectHide()"
     onkeydown="return searchBox.OnSearchSelectKey(event)">
</div>

<!-- iframe showing the search results (closed by default) -->
<div id="MSearchResultsWindow">
<iframe src="javascript:void(0)" frameborder="0" 
        name="MSearchResults" id="MSearchResults">
</iframe>
</div>

<div class="header">
  <div class="headertitle">
<div class="title">Sparse matrix manipulations<div class="ingroups"><a class="el" href="group__Sparse__chapter.html">Sparse linear algebra</a></div></div>  </div>
</div><!--header-->
<div class="contents">
<p>Manipulating and solving sparse problems involves various modules which are summarized below:</p>
<table class="manual">
<tr>
<th>Module</th><th>Header file</th><th>Contents </th></tr>
<tr>
<td><a class="el" href="group__SparseCore__Module.html">SparseCore </a></td><td><div class="fragment"><div class="line"><span class="preprocessor">#include &lt;Eigen/SparseCore&gt;</span></div>
</div><!-- fragment --></td><td><a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation.">SparseMatrix</a> and <a class="el" href="classEigen_1_1SparseVector.html" title="a sparse vector class">SparseVector</a> classes, matrix assembly, basic sparse linear algebra (including sparse triangular solvers) </td></tr>
<tr>
<td><a class="el" href="group__SparseCholesky__Module.html">SparseCholesky </a></td><td><div class="fragment"><div class="line"><span class="preprocessor">#include &lt;Eigen/SparseCholesky&gt;</span></div>
</div><!-- fragment --></td><td>Direct sparse <a class="el" href="classEigen_1_1LLT.html" title="Standard Cholesky decomposition (LL^T) of a matrix and associated features.">LLT</a> and <a class="el" href="classEigen_1_1LDLT.html" title="Robust Cholesky decomposition of a matrix with pivoting.">LDLT</a> Cholesky factorization to solve sparse self-adjoint positive definite problems </td></tr>
<tr>
<td><a class="el" href="group__SparseLU__Module.html">SparseLU </a></td><td><div class="fragment"><div class="line"><span class="preprocessor">#include&lt;Eigen/SparseLU&gt;</span> </div>
</div><!-- fragment --> </td><td>Sparse LU factorization to solve general square sparse systems </td></tr>
<tr>
<td><a class="el" href="group__SparseQR__Module.html">SparseQR </a></td><td><div class="fragment"><div class="line"><span class="preprocessor">#include&lt;Eigen/SparseQR&gt;</span></div>
</div><!-- fragment --> </td><td>Sparse QR factorization for solving sparse linear least-squares problems </td></tr>
<tr>
<td><a class="el" href="group__IterativeLinearSolvers__Module.html">IterativeLinearSolvers </a></td><td><div class="fragment"><div class="line"><span class="preprocessor">#include &lt;Eigen/IterativeLinearSolvers&gt;</span></div>
</div><!-- fragment --></td><td>Iterative solvers to solve large general linear square problems (including self-adjoint positive definite problems) </td></tr>
<tr>
<td><a class="el" href="group__Sparse__Module.html">Sparse </a></td><td><div class="fragment"><div class="line"><span class="preprocessor">#include &lt;Eigen/Sparse&gt;</span></div>
</div><!-- fragment --></td><td>Includes all the above modules </td></tr>
</table>
<h1><a class="anchor" id="TutorialSparseIntro"></a>
Sparse matrix format</h1>
<p>In many applications (e.g., finite element methods) it is common to deal with very large matrices where only a few coefficients are different from zero. In such cases, memory consumption can be reduced and performance increased by using a specialized representation storing only the nonzero coefficients. Such a matrix is called a sparse matrix.</p>
<p><b>The</b> <b>SparseMatrix</b> <b>class</b> </p>
<p>The class <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation.">SparseMatrix</a> is the main sparse matrix representation of <a class="el" href="namespaceEigen.html" title="Namespace containing all symbols from the Eigen library.">Eigen</a>'s sparse module; it offers high performance and low memory usage. It implements a more versatile variant of the widely-used Compressed Column (or Row) Storage scheme. It consists of four compact arrays:</p><ul>
<li><code>Values:</code> stores the coefficient values of the non-zeros.</li>
<li><code>InnerIndices:</code> stores the row (resp. column) indices of the non-zeros.</li>
<li><code>OuterStarts:</code> stores for each column (resp. row) the index of the first non-zero in the previous two arrays.</li>
<li><code>InnerNNZs:</code> stores the number of non-zeros of each column (resp. row). The word <code>inner</code> refers to an <em>inner</em> <em>vector</em> that is a column for a column-major matrix, or a row for a row-major matrix. The word <code>outer</code> refers to the other direction.</li>
</ul>
<p>This storage scheme is better explained on an example. The following matrix </p><table class="manual">
<tr>
<td>0</td><td>3</td><td>0</td><td>0</td><td>0 </td></tr>
<tr>
<td>22</td><td>0</td><td>0</td><td>0</td><td>17 </td></tr>
<tr>
<td>7</td><td>5</td><td>0</td><td>1</td><td>0 </td></tr>
<tr>
<td>0</td><td>0</td><td>0</td><td>0</td><td>0 </td></tr>
<tr>
<td>0</td><td>0</td><td>14</td><td>0</td><td>8 </td></tr>
</table>
<p>and one of its possible sparse, <b>column</b> <b>major</b> representation: </p><table class="manual">
<tr>
<td>Values: </td><td>22</td><td>7</td><td>_</td><td>3</td><td>5</td><td>14</td><td>_</td><td>_</td><td>1</td><td>_</td><td>17</td><td>8 </td></tr>
<tr>
<td>InnerIndices: </td><td>1</td><td>2</td><td>_</td><td>0</td><td>2</td><td>4</td><td>_</td><td>_</td><td>2</td><td>_</td><td>1</td><td>4 </td></tr>
</table>
<table class="manual">
<tr>
<td>OuterStarts:</td><td>0</td><td>3</td><td>5</td><td>8</td><td>10</td><td><em>12</em>  </td></tr>
<tr>
<td>InnerNNZs: </td><td>2</td><td>2</td><td>1</td><td>1</td><td>2</td><td></td></tr>
</table>
<p>Currently the elements of a given inner vector are guaranteed to be always sorted by increasing inner indices. The <code>"_"</code> indicates available free space to quickly insert new elements. Assuming no reallocation is needed, the insertion of a random element is therefore in <code>O(nnz_j)</code> where <code>nnz_j</code> is the number of nonzeros of the respective inner vector. On the other hand, inserting elements with increasing inner indices in a given inner vector is much more efficient since this only requires to increase the respective <code>InnerNNZs</code> entry that is a <code>O(1)</code> operation.</p>
<p>The case where no empty space is available is a special case, and is referred as the <em>compressed</em> mode. It corresponds to the widely used Compressed Column (or Row) Storage schemes (CCS or CRS). Any <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation.">SparseMatrix</a> can be turned to this form by calling the <a class="el" href="classEigen_1_1SparseMatrix.html#a6400e6e733a7131f9b776d6386d2ed89">SparseMatrix::makeCompressed()</a> function. In this case, one can remark that the <code>InnerNNZs</code> array is redundant with <code>OuterStarts</code> because we have the equality: <code>InnerNNZs[j] == OuterStarts[j+1] - OuterStarts[j]</code>. Therefore, in practice a call to <a class="el" href="classEigen_1_1SparseMatrix.html#a6400e6e733a7131f9b776d6386d2ed89">SparseMatrix::makeCompressed()</a> frees this buffer.</p>
<p>It is worth noting that most of our wrappers to external libraries requires compressed matrices as inputs.</p>
<p>The results of Eigen's operations always produces <b>compressed</b> sparse matrices. On the other hand, the insertion of a new element into a <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation.">SparseMatrix</a> converts this later to the <b>uncompressed</b> mode.</p>
<p>Here is the previous matrix represented in compressed mode: </p><table class="manual">
<tr>
<td>Values: </td><td>22</td><td>7</td><td>3</td><td>5</td><td>14</td><td>1</td><td>17</td><td>8 </td></tr>
<tr>
<td>InnerIndices: </td><td>1</td><td>2</td><td>0</td><td>2</td><td>4</td><td>2</td><td>1</td><td>4 </td></tr>
</table>
<table class="manual">
<tr>
<td>OuterStarts:</td><td>0</td><td>2</td><td>4</td><td>5</td><td>6</td><td><em>8</em>  </td></tr>
</table>
<p>A <a class="el" href="classEigen_1_1SparseVector.html" title="a sparse vector class">SparseVector</a> is a special case of a <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation.">SparseMatrix</a> where only the <code>Values</code> and <code>InnerIndices</code> arrays are stored. There is no notion of compressed/uncompressed mode for a <a class="el" href="classEigen_1_1SparseVector.html" title="a sparse vector class">SparseVector</a>.</p>
<h1><a class="anchor" id="TutorialSparseExample"></a>
First example</h1>
<p>Before describing each individual class, let's start with the following typical example: solving the Laplace equation \( \Delta u = 0 \) on a regular 2D grid using a finite difference scheme and Dirichlet boundary conditions. Such problem can be mathematically expressed as a linear problem of the form \( Ax=b \) where \( x \) is the vector of <code>m</code> unknowns (in our case, the values of the pixels), \( b \) is the right hand side vector resulting from the boundary conditions, and \( A \) is an \( m \times m \) matrix containing only a few non-zero elements resulting from the discretization of the Laplacian operator.</p>
<table class="manual">
<tr>
<td><div class="fragment"><div class="line"><span class="preprocessor">#include &lt;Eigen/Sparse&gt;</span></div>
<div class="line"><span class="preprocessor">#include &lt;vector&gt;</span></div>
<div class="line"><span class="preprocessor">#include &lt;iostream&gt;</span></div>
<div class="line"> </div>
<div class="line"><span class="keyword">typedef</span> <a class="code" href="classEigen_1_1SparseMatrix.html">Eigen::SparseMatrix&lt;double&gt;</a> SpMat; <span class="comment">// declares a column-major sparse matrix type of double</span></div>
<div class="line"><span class="keyword">typedef</span> <a class="code" href="classEigen_1_1Triplet.html">Eigen::Triplet&lt;double&gt;</a> T;</div>
<div class="line"> </div>
<div class="line"><span class="keywordtype">void</span> buildProblem(std::vector&lt;T&gt;&amp; coefficients, <a class="code" href="classEigen_1_1Matrix.html">Eigen::VectorXd</a>&amp; b, <span class="keywordtype">int</span> n);</div>
<div class="line"><span class="keywordtype">void</span> saveAsBitmap(<span class="keyword">const</span> <a class="code" href="classEigen_1_1Matrix.html">Eigen::VectorXd</a>&amp; x, <span class="keywordtype">int</span> n, <span class="keyword">const</span> <span class="keywordtype">char</span>* filename);</div>
<div class="line"> </div>
<div class="line"><span class="keywordtype">int</span> main(<span class="keywordtype">int</span> argc, <span class="keywordtype">char</span>** argv)</div>
<div class="line">{</div>
<div class="line">  <span class="keywordflow">if</span>(argc!=2) {</div>
<div class="line">    std::cerr &lt;&lt; <span class="stringliteral">&quot;Error: expected one and only one argument.\n&quot;</span>;</div>
<div class="line">    <span class="keywordflow">return</span> -1;</div>
<div class="line">  }</div>
<div class="line">  </div>
<div class="line">  <span class="keywordtype">int</span> n = 300;  <span class="comment">// size of the image</span></div>
<div class="line">  <span class="keywordtype">int</span> m = n*n;  <span class="comment">// number of unknowns (=number of pixels)</span></div>
<div class="line"> </div>
<div class="line">  <span class="comment">// Assembly:</span></div>
<div class="line">  std::vector&lt;T&gt; coefficients;            <span class="comment">// list of non-zeros coefficients</span></div>
<div class="line">  <a class="code" href="classEigen_1_1Matrix.html">Eigen::VectorXd</a> b(m);                   <span class="comment">// the right hand side-vector resulting from the constraints</span></div>
<div class="line">  buildProblem(coefficients, b, n);</div>
<div class="line"> </div>
<div class="line">  SpMat A(m,m);</div>
<div class="line">  A.setFromTriplets(coefficients.begin(), coefficients.end());</div>
<div class="line"> </div>
<div class="line">  <span class="comment">// Solving:</span></div>
<div class="line">  <a class="code" href="classEigen_1_1SimplicialCholesky.html">Eigen::SimplicialCholesky&lt;SpMat&gt;</a> chol(A);  <span class="comment">// performs a Cholesky factorization of A</span></div>
<div class="line">  <a class="code" href="classEigen_1_1Matrix.html">Eigen::VectorXd</a> x = chol.solve(b);         <span class="comment">// use the factorization to solve for the given right hand side</span></div>
<div class="line"> </div>
<div class="line">  <span class="comment">// Export the result to a file:</span></div>
<div class="line">  saveAsBitmap(x, n, argv[1]);</div>
<div class="line"> </div>
<div class="line">  <span class="keywordflow">return</span> 0;</div>
<div class="line">}</div>
<div class="line"> </div>
<div class="ttc" id="aclassEigen_1_1Matrix_html"><div class="ttname"><a href="classEigen_1_1Matrix.html">Eigen::Matrix</a></div><div class="ttdoc">The matrix class, also used for vectors and row-vectors.</div><div class="ttdef"><b>Definition:</b> Matrix.h:182</div></div>
<div class="ttc" id="aclassEigen_1_1SimplicialCholesky_html"><div class="ttname"><a href="classEigen_1_1SimplicialCholesky.html">Eigen::SimplicialCholesky</a></div><div class="ttdef"><b>Definition:</b> SimplicialCholesky.h:514</div></div>
<div class="ttc" id="aclassEigen_1_1SparseMatrix_html"><div class="ttname"><a href="classEigen_1_1SparseMatrix.html">Eigen::SparseMatrix</a></div><div class="ttdoc">A versatible sparse matrix representation.</div><div class="ttdef"><b>Definition:</b> SparseMatrix.h:100</div></div>
<div class="ttc" id="aclassEigen_1_1Triplet_html"><div class="ttname"><a href="classEigen_1_1Triplet.html">Eigen::Triplet</a></div><div class="ttdoc">A small structure to hold a non zero as a triplet (i,j,value).</div><div class="ttdef"><b>Definition:</b> SparseUtil.h:163</div></div>
</div><!-- fragment -->  </td><td><div class="image">
<img src="Tutorial_sparse_example.jpeg" alt=""/>
</div>
 </td></tr>
</table>
<p>In this example, we start by defining a column-major sparse matrix type of double <code>SparseMatrix&lt;double&gt;</code>, and a triplet list of the same scalar type <code>Triplet&lt;double&gt;</code>. A triplet is a simple object representing a non-zero entry as the triplet: <code>row</code> index, <code>column</code> index, <code>value</code>.</p>
<p>In the main function, we declare a list <code>coefficients</code> of triplets (as a std vector) and the right hand side vector \( b \) which are filled by the <em>buildProblem</em> function. The raw and flat list of non-zero entries is then converted to a true <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation.">SparseMatrix</a> object <code>A</code>. Note that the elements of the list do not have to be sorted, and possible duplicate entries will be summed up.</p>
<p>The last step consists of effectively solving the assembled problem. Since the resulting matrix <code>A</code> is symmetric by construction, we can perform a direct Cholesky factorization via the <a class="el" href="classEigen_1_1SimplicialLDLT.html" title="A direct sparse LDLT Cholesky factorizations without square root.">SimplicialLDLT</a> class which behaves like its <a class="el" href="classEigen_1_1LDLT.html" title="Robust Cholesky decomposition of a matrix with pivoting.">LDLT</a> counterpart for dense objects.</p>
<p>The resulting vector <code>x</code> contains the pixel values as a 1D array which is saved to a jpeg file shown on the right of the code above.</p>
<p>Describing the <em>buildProblem</em> and <em>save</em> functions is out of the scope of this tutorial. They are given <a class="el" href="TutorialSparse_example_details.html">here</a> for the curious and reproducibility purpose.</p>
<h1><a class="anchor" id="TutorialSparseSparseMatrix"></a>
The SparseMatrix class</h1>
<p><b>Matrix</b> <b>and</b> <b>vector</b> <b>properties</b> <br  />
 The <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation.">SparseMatrix</a> and <a class="el" href="classEigen_1_1SparseVector.html" title="a sparse vector class">SparseVector</a> classes take three template arguments: the scalar type (e.g., double) the storage order (ColMajor or RowMajor, the default is ColMajor) the inner index type (default is <code>int</code>).</p>
<p>As for dense <a class="el" href="classEigen_1_1Matrix.html" title="The matrix class, also used for vectors and row-vectors.">Matrix</a> objects, constructors takes the size of the object. Here are some examples:</p>
<div class="fragment"><div class="line">SparseMatrix&lt;std::complex&lt;float&gt; &gt; mat(1000,2000);         <span class="comment">// declares a 1000x2000 column-major compressed sparse matrix of complex&lt;float&gt;</span></div>
<div class="line">SparseMatrix&lt;double,RowMajor&gt; mat(1000,2000);              <span class="comment">// declares a 1000x2000 row-major compressed sparse matrix of double</span></div>
<div class="line">SparseVector&lt;std::complex&lt;float&gt; &gt; vec(1000);              <span class="comment">// declares a column sparse vector of complex&lt;float&gt; of size 1000</span></div>
<div class="line">SparseVector&lt;double,RowMajor&gt; vec(1000);                   <span class="comment">// declares a row sparse vector of double of size 1000</span></div>
</div><!-- fragment --><p>In the rest of the tutorial, <code>mat</code> and <code>vec</code> represent any sparse-matrix and sparse-vector objects, respectively.</p>
<p>The dimensions of a matrix can be queried using the following functions: </p><table class="manual">
<tr>
<td>Standard <br  />
 dimensions</td><td><div class="fragment"><div class="line">mat.rows()</div>
<div class="line">mat.cols()</div>
</div><!-- fragment --> </td><td><div class="fragment"><div class="line">vec.size() </div>
</div><!-- fragment -->  </td></tr>
<tr>
<td>Sizes along the <br  />
 inner/outer dimensions</td><td><div class="fragment"><div class="line">mat.innerSize()</div>
<div class="line">mat.outerSize()</div>
</div><!-- fragment --> </td><td></td></tr>
<tr>
<td>Number of non <br  />
 zero coefficients</td><td><div class="fragment"><div class="line">mat.nonZeros() </div>
</div><!-- fragment --> </td><td><div class="fragment"><div class="line">vec.nonZeros() </div>
</div><!-- fragment --> </td></tr>
</table>
<p><b>Iterating</b> <b>over</b> <b>the</b> <b>nonzero</b> <b>coefficients</b> <br  />
 Random access to the elements of a sparse object can be done through the <code>coeffRef(i,j)</code> function. However, this function involves a quite expensive binary search. In most cases, one only wants to iterate over the non-zeros elements. This is achieved by a standard loop over the outer dimension, and then by iterating over the non-zeros of the current inner vector via an InnerIterator. Thus, the non-zero entries have to be visited in the same order than the storage order. Here is an example: </p><table class="manual">
<tr>
<td><div class="fragment"><div class="line">SparseMatrix&lt;double&gt; mat(rows,cols);</div>
<div class="line"><span class="keywordflow">for</span> (<span class="keywordtype">int</span> k=0; k&lt;mat.outerSize(); ++k)</div>
<div class="line">  <span class="keywordflow">for</span> (SparseMatrix&lt;double&gt;::InnerIterator it(mat,k); it; ++it)</div>
<div class="line">  {</div>
<div class="line">    it.value();</div>
<div class="line">    it.row();   <span class="comment">// row index</span></div>
<div class="line">    it.col();   <span class="comment">// col index (here it is equal to k)</span></div>
<div class="line">    it.index(); <span class="comment">// inner index, here it is equal to it.row()</span></div>
<div class="line">  }</div>
</div><!-- fragment --> </td><td><div class="fragment"><div class="line">SparseVector&lt;double&gt; vec(size);</div>
<div class="line"><span class="keywordflow">for</span> (SparseVector&lt;double&gt;::InnerIterator it(vec); it; ++it)</div>
<div class="line">{</div>
<div class="line">  it.value(); <span class="comment">// == vec[ it.index() ]</span></div>
<div class="line">  it.index();</div>
<div class="line">}</div>
</div><!-- fragment -->  </td></tr>
</table>
<p>For a writable expression, the referenced value can be modified using the valueRef() function. If the type of the sparse matrix or vector depends on a template parameter, then the <code>typename</code> keyword is required to indicate that <code>InnerIterator</code> denotes a type; see <a class="el" href="TopicTemplateKeyword.html">The template and typename keywords in C++</a> for details.</p>
<h1><a class="anchor" id="TutorialSparseFilling"></a>
Filling a sparse matrix</h1>
<p>Because of the special storage scheme of a <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation.">SparseMatrix</a>, special care has to be taken when adding new nonzero entries. For instance, the cost of a single purely random insertion into a <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation.">SparseMatrix</a> is <code>O(nnz)</code>, where <code>nnz</code> is the current number of non-zero coefficients.</p>
<p>The simplest way to create a sparse matrix while guaranteeing good performance is thus to first build a list of so-called <em>triplets</em>, and then convert it to a <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation.">SparseMatrix</a>.</p>
<p>Here is a typical usage example: </p><div class="fragment"><div class="line"><span class="keyword">typedef</span> <a class="code" href="classEigen_1_1Triplet.html">Eigen::Triplet&lt;double&gt;</a> T;</div>
<div class="line">std::vector&lt;T&gt; tripletList;</div>
<div class="line">tripletList.reserve(estimation_of_entries);</div>
<div class="line"><span class="keywordflow">for</span>(...)</div>
<div class="line">{</div>
<div class="line">  <span class="comment">// ...</span></div>
<div class="line">  tripletList.push_back(T(i,j,v_ij));</div>
<div class="line">}</div>
<div class="line">SparseMatrixType mat(rows,cols);</div>
<div class="line">mat.setFromTriplets(tripletList.begin(), tripletList.end());</div>
<div class="line"><span class="comment">// mat is ready to go!</span></div>
</div><!-- fragment --><p> The <code>std::vector</code> of triplets might contain the elements in arbitrary order, and might even contain duplicated elements that will be summed up by setFromTriplets(). See the <a class="el" href="classEigen_1_1SparseMatrix.html#a8f09e3597f37aa8861599260af6a53e0">SparseMatrix::setFromTriplets()</a> function and class <a class="el" href="classEigen_1_1Triplet.html" title="A small structure to hold a non zero as a triplet (i,j,value).">Triplet</a> for more details.</p>
<p>In some cases, however, slightly higher performance, and lower memory consumption can be reached by directly inserting the non-zeros into the destination matrix. A typical scenario of this approach is illustrated below: </p><div class="fragment"><div class="line">1: SparseMatrix&lt;double&gt; mat(rows,cols);         <span class="comment">// default is column major</span></div>
<div class="line">2: mat.reserve(<a class="code" href="classEigen_1_1DenseBase.html#a68a7ece6c5629d1e9447a321fcb14ccd">VectorXi::Constant</a>(cols,6));</div>
<div class="line">3: <span class="keywordflow">for each</span> i,j such that v_ij != 0</div>
<div class="line">4:   mat.insert(i,j) = v_ij;                    <span class="comment">// alternative: mat.coeffRef(i,j) += v_ij;</span></div>
<div class="line">5: mat.makeCompressed();                        <span class="comment">// optional</span></div>
<div class="ttc" id="aclassEigen_1_1DenseBase_html_a68a7ece6c5629d1e9447a321fcb14ccd"><div class="ttname"><a href="classEigen_1_1DenseBase.html#a68a7ece6c5629d1e9447a321fcb14ccd">Eigen::DenseBase::Constant</a></div><div class="ttdeci">static const ConstantReturnType Constant(Index rows, Index cols, const Scalar &amp;value)</div><div class="ttdef"><b>Definition:</b> CwiseNullaryOp.h:191</div></div>
</div><!-- fragment --><ul>
<li>The key ingredient here is the line 2 where we reserve room for 6 non-zeros per column. In many cases, the number of non-zeros per column or row can easily be known in advance. If it varies significantly for each inner vector, then it is possible to specify a reserve size for each inner vector by providing a vector object with an <code>operator[](int j)</code> returning the reserve size of the <code>j-th</code> inner vector (e.g., via a <code>VectorXi</code> or <code>std::vector&lt;int&gt;</code>). If only a rought estimate of the number of nonzeros per inner-vector can be obtained, it is highly recommended to overestimate it rather than the opposite. If this line is omitted, then the first insertion of a new element will reserve room for 2 elements per inner vector.</li>
<li>The line 4 performs a sorted insertion. In this example, the ideal case is when the <code>j-th</code> column is not full and contains non-zeros whose inner-indices are smaller than <code>i</code>. In this case, this operation boils down to trivial O(1) operation.</li>
<li>When calling <code>insert(i,j)</code> the element <code>i</code>, <code>j</code> must not already exists, otherwise use the <code>coeffRef(i,j)</code> method that will allow to, e.g., accumulate values. This method first performs a binary search and finally calls <code>insert(i,j)</code> if the element does not already exist. It is more flexible than <code>insert()</code> but also more costly.</li>
<li>The line 5 suppresses the remaining empty space and transforms the matrix into a compressed column storage.</li>
</ul>
<h1><a class="anchor" id="TutorialSparseFeatureSet"></a>
Supported operators and functions</h1>
<p>Because of their special storage format, sparse matrices cannot offer the same level of flexibility than dense matrices. In <a class="el" href="namespaceEigen.html" title="Namespace containing all symbols from the Eigen library.">Eigen</a>'s sparse module we chose to expose only the subset of the dense matrix API which can be efficiently implemented. In the following <em>sm</em> denotes a sparse matrix, <em>sv</em> a sparse vector, <em>dm</em> a dense matrix, and <em>dv</em> a dense vector.</p>
<h2><a class="anchor" id="TutorialSparse_BasicOps"></a>
Basic operations</h2>
<p>Sparse expressions support most of the unary and binary coefficient wise operations: </p><div class="fragment"><div class="line">sm1.real()   sm1.imag()   -sm1                    0.5*sm1</div>
<div class="line">sm1+sm2      sm1-sm2      sm1.cwiseProduct(sm2)</div>
</div><!-- fragment --><p> However, <b>a strong restriction is that the storage orders must match</b>. For instance, in the following example: </p><div class="fragment"><div class="line">sm4 = sm1 + sm2 + sm3;</div>
</div><!-- fragment --><p> sm1, sm2, and sm3 must all be row-major or all column-major. On the other hand, there is no restriction on the target matrix sm4. For instance, this means that for computing \( A^T + A \), the matrix \( A^T \) must be evaluated into a temporary matrix of compatible storage order: </p><div class="fragment"><div class="line">SparseMatrix&lt;double&gt; A, B;</div>
<div class="line">B = SparseMatrix&lt;double&gt;(A.transpose()) + A;</div>
</div><!-- fragment --><p>Binary coefficient wise operators can also mix sparse and dense expressions: </p><div class="fragment"><div class="line">sm2 = sm1.cwiseProduct(dm1);</div>
<div class="line">dm2 = sm1 + dm1;</div>
<div class="line">dm2 = dm1 - sm1;</div>
</div><!-- fragment --><p> Performance-wise, the adding/subtracting sparse and dense matrices is better performed in two steps. For instance, instead of doing <code>dm2 = sm1 + dm1</code>, better write: </p><div class="fragment"><div class="line">dm2 = dm1;</div>
<div class="line">dm2 += sm1;</div>
</div><!-- fragment --><p> This version has the advantage to fully exploit the higher performance of dense storage (no indirection, SIMD, etc.), and to pay the cost of slow sparse evaluation on the few non-zeros of the sparse matrix only.</p>
<p>Sparse expressions also support transposition: </p><div class="fragment"><div class="line">sm1 = sm2.transpose();</div>
<div class="line">sm1 = sm2.adjoint();</div>
</div><!-- fragment --><p> However, there is no <code>transposeInPlace()</code> method.</p>
<h2><a class="anchor" id="TutorialSparse_Products"></a>
Matrix products</h2>
<p>Eigen supports various kind of sparse matrix products which are summarize below:</p><ul>
<li><b>sparse-dense</b>: <div class="fragment"><div class="line">dv2 = sm1 * dv1;</div>
<div class="line">dm2 = dm1 * sm1.adjoint();</div>
<div class="line">dm2 = 2. * sm1 * dm1;</div>
</div><!-- fragment --></li>
<li><b>symmetric</b> <b>sparse-dense</b>. The product of a sparse symmetric matrix with a dense matrix (or vector) can also be optimized by specifying the symmetry with <code>selfadjointView()</code>: <div class="fragment"><div class="line">dm2 = sm1.selfadjointView&lt;&gt;() * dm1;          <span class="comment">// if all coefficients of sm1 are stored</span></div>
<div class="line">dm2 = sm1.selfadjointView&lt;<a class="code" href="group__enums.html#gga39e3366ff5554d731e7dc8bb642f83cdafca2ccebb604f171656deb53e8c083c1">Upper</a>&gt;() * dm1;     <span class="comment">// if only the upper part of sm1 is stored</span></div>
<div class="line">dm2 = sm1.selfadjointView&lt;<a class="code" href="group__enums.html#gga39e3366ff5554d731e7dc8bb642f83cdaf581029282d421eee5aae14238c6f749">Lower</a>&gt;() * dm1;     <span class="comment">// if only the lower part of sm1 is stored</span></div>
<div class="ttc" id="agroup__enums_html_gga39e3366ff5554d731e7dc8bb642f83cdaf581029282d421eee5aae14238c6f749"><div class="ttname"><a href="group__enums.html#gga39e3366ff5554d731e7dc8bb642f83cdaf581029282d421eee5aae14238c6f749">Eigen::Lower</a></div><div class="ttdeci">@ Lower</div><div class="ttdef"><b>Definition:</b> Constants.h:211</div></div>
<div class="ttc" id="agroup__enums_html_gga39e3366ff5554d731e7dc8bb642f83cdafca2ccebb604f171656deb53e8c083c1"><div class="ttname"><a href="group__enums.html#gga39e3366ff5554d731e7dc8bb642f83cdafca2ccebb604f171656deb53e8c083c1">Eigen::Upper</a></div><div class="ttdeci">@ Upper</div><div class="ttdef"><b>Definition:</b> Constants.h:213</div></div>
</div><!-- fragment --></li>
<li><b>sparse-sparse</b>. For sparse-sparse products, two different algorithms are available. The default one is conservative and preserve the explicit zeros that might appear: <div class="fragment"><div class="line">sm3 = sm1 * sm2;</div>
<div class="line">sm3 = 4 * sm1.adjoint() * sm2;</div>
</div><!-- fragment --> The second algorithm prunes on the fly the explicit zeros, or the values smaller than a given threshold. It is enabled and controlled through the <code>prune()</code> functions: <div class="fragment"><div class="line">sm3 = (sm1 * sm2).pruned();                  <span class="comment">// removes numerical zeros</span></div>
<div class="line">sm3 = (sm1 * sm2).pruned(ref);               <span class="comment">// removes elements much smaller than ref</span></div>
<div class="line">sm3 = (sm1 * sm2).pruned(ref,epsilon);       <span class="comment">// removes elements smaller than ref*epsilon</span></div>
</div><!-- fragment --></li>
<li><b>permutations</b>. Finally, permutations can be applied to sparse matrices too: <div class="fragment"><div class="line">PermutationMatrix&lt;Dynamic,Dynamic&gt; P = ...;</div>
<div class="line">sm2 = P * sm1;</div>
<div class="line">sm2 = sm1 * P.inverse();</div>
<div class="line">sm2 = sm1.transpose() * P;</div>
</div><!-- fragment --></li>
</ul>
<h2><a class="anchor" id="TutorialSparse_SubMatrices"></a>
Block operations</h2>
<p>Regarding read-access, sparse matrices expose the same API than for dense matrices to access to sub-matrices such as blocks, columns, and rows. See <a class="el" href="group__TutorialBlockOperations.html">Block operations</a> for a detailed introduction. However, for performance reasons, writing to a sub-sparse-matrix is much more limited, and currently only contiguous sets of columns (resp. rows) of a column-major (resp. row-major) <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation.">SparseMatrix</a> are writable. Moreover, this information has to be known at compile-time, leaving out methods such as <code>block(...)</code> and <code>corner*(...)</code>. The available API for write-access to a <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation.">SparseMatrix</a> are summarized below: </p><div class="fragment"><div class="line">SparseMatrix&lt;double,ColMajor&gt; sm1;</div>
<div class="line">sm1.col(j) = ...;</div>
<div class="line">sm1.leftCols(ncols) = ...;</div>
<div class="line">sm1.middleCols(j,ncols) = ...;</div>
<div class="line">sm1.rightCols(ncols) = ...;</div>
<div class="line"> </div>
<div class="line">SparseMatrix&lt;double,RowMajor&gt; sm2;</div>
<div class="line">sm2.row(i) = ...;</div>
<div class="line">sm2.topRows(nrows) = ...;</div>
<div class="line">sm2.middleRows(i,nrows) = ...;</div>
<div class="line">sm2.bottomRows(nrows) = ...;</div>
</div><!-- fragment --><p>In addition, sparse matrices expose the <code>SparseMatrixBase::innerVector()</code> and <code>SparseMatrixBase::innerVectors()</code> methods, which are aliases to the <code>col</code>/<code>middleCols</code> methods for a column-major storage, and to the <code>row</code>/<code>middleRows</code> methods for a row-major storage.</p>
<h2><a class="anchor" id="TutorialSparse_TriangularSelfadjoint"></a>
Triangular and selfadjoint views</h2>
<p>Just as with dense matrices, the <code>triangularView()</code> function can be used to address a triangular part of the matrix, and perform triangular solves with a dense right hand side: </p><div class="fragment"><div class="line">dm2 = sm1.triangularView&lt;<a class="code" href="group__enums.html#gga39e3366ff5554d731e7dc8bb642f83cdaf581029282d421eee5aae14238c6f749">Lower</a>&gt;(dm1);</div>
<div class="line">dv2 = sm1.transpose().triangularView&lt;<a class="code" href="group__enums.html#gga39e3366ff5554d731e7dc8bb642f83cdafca2ccebb604f171656deb53e8c083c1">Upper</a>&gt;(dv1);</div>
</div><!-- fragment --><p>The <code>selfadjointView()</code> function permits various operations:</p><ul>
<li>optimized sparse-dense matrix products: <div class="fragment"><div class="line">dm2 = sm1.selfadjointView&lt;&gt;() * dm1;          <span class="comment">// if all coefficients of sm1 are stored</span></div>
<div class="line">dm2 = sm1.selfadjointView&lt;<a class="code" href="group__enums.html#gga39e3366ff5554d731e7dc8bb642f83cdafca2ccebb604f171656deb53e8c083c1">Upper</a>&gt;() * dm1;     <span class="comment">// if only the upper part of sm1 is stored</span></div>
<div class="line">dm2 = sm1.selfadjointView&lt;<a class="code" href="group__enums.html#gga39e3366ff5554d731e7dc8bb642f83cdaf581029282d421eee5aae14238c6f749">Lower</a>&gt;() * dm1;     <span class="comment">// if only the lower part of sm1 is stored</span></div>
</div><!-- fragment --></li>
<li>copy of triangular parts: <div class="fragment"><div class="line">sm2 = sm1.selfadjointView&lt;<a class="code" href="group__enums.html#gga39e3366ff5554d731e7dc8bb642f83cdafca2ccebb604f171656deb53e8c083c1">Upper</a>&gt;();                               <span class="comment">// makes a full selfadjoint matrix from the upper triangular part</span></div>
<div class="line">sm2.selfadjointView&lt;<a class="code" href="group__enums.html#gga39e3366ff5554d731e7dc8bb642f83cdaf581029282d421eee5aae14238c6f749">Lower</a>&gt;() = sm1.selfadjointView&lt;<a class="code" href="group__enums.html#gga39e3366ff5554d731e7dc8bb642f83cdafca2ccebb604f171656deb53e8c083c1">Upper</a>&gt;();      <span class="comment">// copies the upper triangular part to the lower triangular part</span></div>
</div><!-- fragment --></li>
<li>application of symmetric permutations: <div class="fragment"><div class="line">PermutationMatrix&lt;Dynamic,Dynamic&gt; P = ...;</div>
<div class="line">sm2 = A.selfadjointView&lt;<a class="code" href="group__enums.html#gga39e3366ff5554d731e7dc8bb642f83cdafca2ccebb604f171656deb53e8c083c1">Upper</a>&gt;().twistedBy(P);                                <span class="comment">// compute P S P&#39; from the upper triangular part of A, and make it a full matrix</span></div>
<div class="line">sm2.selfadjointView&lt;<a class="code" href="group__enums.html#gga39e3366ff5554d731e7dc8bb642f83cdaf581029282d421eee5aae14238c6f749">Lower</a>&gt;() = A.selfadjointView&lt;<a class="code" href="group__enums.html#gga39e3366ff5554d731e7dc8bb642f83cdaf581029282d421eee5aae14238c6f749">Lower</a>&gt;().twistedBy(P);       <span class="comment">// compute P S P&#39; from the lower triangular part of A, and then only compute the lower part</span></div>
</div><!-- fragment --></li>
</ul>
<p>Please, refer to the <a class="el" href="group__SparseQuickRefPage.html">Quick Reference </a> guide for the list of supported operations. The list of linear solvers available is <a class="el" href="group__TopicSparseSystems.html">here. </a> </p>
</div><!-- contents -->
</div><!-- doc-content -->
<!-- start footer part -->
<div id="nav-path" class="navpath"><!-- id is needed for treeview function! -->
  <ul>
    <li class="footer">Generated on Thu Apr 21 2022 13:07:55 for Eigen by
    <a href="http://www.doxygen.org/index.html">
    <img class="footer" src="doxygen.png" alt="doxygen"/></a> 1.9.1 </li>
  </ul>
</div>
</body>
</html>
